CONTRABASS: exploiting flux constraints in genome-scale models for the detection of vulnerabilities

Abstract Motivation Despite the fact that antimicrobial resistance is an increasing health concern, the pace of production of new drugs is slow due to the high cost and uncertain success of the process. The development of high-throughput technologies has allowed the integration of biological data into detailed genome-scale models of multiple organisms. Such models can be exploited by means of computational methods to identify system vulnerabilities such as chokepoint reactions and essential reactions. These vulnerabilities are appealing drug targets that can lead to novel drug developments. However, the current approach to compute these vulnerabilities is only based on topological data and ignores the dynamic information of the model. This can lead to misidentified drug targets. Results This work computes flux constraints that are consistent with a certain growth rate of the modelled organism, and integrates the computed flux constraints into the model to improve the detection of vulnerabilities. By exploiting these flux constraints, we are able to obtain a directionality of the reactions of metabolism consistent with a given growth rate of the model, and consequently, a more realistic detection of vulnerabilities can be performed. Several sets of reactions that are system vulnerabilities are defined and the relationships among them are studied. The approach for the detection of these vulnerabilities has been implemented in the Python tool CONTRABASS. Such tool, for which an online web server has also been implemented, computes flux constraints and generates a report with the detected vulnerabilities. Availability and implementation CONTRABASS is available as an open source Python package at https://github.com/openCONTRABASS/CONTRABASS under GPL-3.0 License. An online web server is available at http://contrabass.unizar.es. Supplementary information A glossary of terms are available at Bioinformatics online.


Introduction
Antimicrobial resistance (AMR) occurs when bacteria, viruses, fungi and parasites evolve to become unresponsive to conventional antibiotics, thus threatening our ability to treat common infections and increasing the risk of spread of severe illnesses. The emergence of multidrug-resistant bacteria (MDR) is particularly alarming, as they can cause infections untreatable with existing antibiotics. According to the World Health Organization, AMR is one of the primary global public health threats of humanity due to its high increasing rate (Leung et al., 2011).
Despite the increasing emergency of AMR, the pipeline of drug development involves a huge cost of time and budget and, in addition, most drug candidates fail the clinical stages. Therefore, novel approaches for the development of drugs that tackle the rising problem of AMR should be established urgently (Jubeh et al., 2020).
The emergence of high-throughput technologies has led to the construction of large datasets with structural and transcriptomic data of different pathogens. The integration of multidimensional data has the potential to offer a more rapid and cost-effective strategy compared to traditional screening methods in the process of drug discovery. In this context, the integration of metabolic functions obtained from genome annotations enabled the reconstruction of metabolic networks of multiple microorganisms. These reconstructions allow the modelling of metabolism, which can provide non-intuitive insights into biological systems that in vivo assays alone cannot provide (Ramos et al., 2018;Richelle et al., 2020).
Metabolism is the set of basic life processes that take place in the cell, and it is the means by which cells can maintain life and grow from their environment. Metabolism can be represented as a metabolic network, which includes all the metabolic reactions that can occur in a cell. A possible strategy for drug development is to find vulnerabilities in the metabolism that can stop the growth and replication of the bacterium, and thus its damaging effects.
As of 2019, genome-scale models (GEMs) of metabolism have been reconstructed for more than 6000 organisms including bacteria, archaea and eukaryotes by either automatic or manual reconstruction (Gu et al., 2019). GEMs have proven useful in a wide range of applications, such as expanding knowledge on microorganisms (Henry et al., 2006;Navid and Almaas, 2012) and microbial communities (Kumar et al., 2019;Zomorrodi and Maranas, 2012), microbial engineering (McAnulty et al., 2012;Roberts et al., 2010) and drug discovery (Ramos et al., 2018). Furthermore, these models have proven to be interesting in areas such as oncology, by studying drug targets in cancer metabolism (Folger et al., 2011), and viral diseases (Bannerman et al., 2021). Drug targeting in pathogens is usually performed by considering essential genes, reactions or metabolites, whose inhibition can effectively kill a pathogen (Gu et al., 2019) or through metabolic network topology analysis (Ramos et al., 2018).
Despite the number of GEMs developed recently, most of these models lack genetic information, which could hamper the model analysis based on gene essentiality and calls for the design of computational methods that exploit as much as possible the available biological information. Here, we focus on the computation of critical reactions in metabolic networks such as essential reactions, dead reactions, and chokepoints (Chukualim et al., 2008;Yeh et al., 2004).
A reaction is essential if it is required by the organism to grow. The deletion of one of these reactions has the potential to cause a malfunction in the metabolism and stop the growth of the organism. As this eventually leads to the death of the organism, essential reactions are widely accepted as drug target candidates (Oyelade et al., 2018).
Dead reactions are those reactions that cannot carry any flux. These reactions will receive special attention since they may reflect an incompleteness of the model or non-preferred pathways of the organism. As it will be shown, these reactions directly affect the number of critical reactions computed for a model.
A chokepoint is a reaction that is the only consumer or the only producer of a given metabolite. Thus, the inhibition of a chokepoint may lead to the unlimited accumulation of potentially toxic metabolites or the lack of production of an essential compound. Chokepoints are, therefore, appealing drug targets in microbial metabolism.
The current methods for the computation of chokepoints only make use of the topological properties of the metabolic network. Although the directionality of the network accounts for the sets of reactants and products of each reaction, it disregards the fact that such sets depend on the fluxes of the reactions.
Reaction fluxes define the direction in which they take place. To establish a particular direction, all reactions are given flux bounds that delimit the range of flux values that a reaction can take. Reactions that are given positive upper flux bounds and negative lower flux bounds are said to be reversible. These flux bounds are used to model reactions that can take place in both directions or simply to model reactions whose direction is unknown (Thiele and Palsson, 2010). For instance, let us say that a chemical reaction r : A $ B is reversible. This means that r can happen in a forward direction (r : A ! B) when a positive flux is assigned, or a backward direction when a negative flux is assigned (r : A B). Notice that this means that metabolites in reversible reactions can act both as reactants and as products. In the reaction r : A ! B, metabolite A is the reactant and metabolite B is the product, however, if the flux of the reaction is negative, then A becomes the product and B the reactant.
In most GEMs, unknown flux bounds are given default flux bounds, e.g. bðrÞ ¼ À1000 mmol gDW À1 h À1 and ubðrÞ ¼ 1000 mmol gDW À1 h À1 , that is, reactions are defined as reversible by default. However, not all these default flux bounds are compatible with a positive growth rate. This work exploits the flux constraints of the model to compute sets of reactants and products that are consistent with a given growth rate, and in turn, to compute growth-dependent critical reactions. Such approach improves the detection of vulnerabilities and is applied to compute chokepoints, essential reactions and dead reactions. The dependence of these sets with respect to growth is also inferred theoretically. To facilitate the computation of vulnerabilities in GEMs by integrating dynamic constraints, a software tool (CONTRABASS) has been developed.
The remainder of the article is organized as follows: Sections 2.1 and 2.2 recall preliminary concepts and definitions such as constraint-based models and essential reactions. Structural definitions of constraint-based models and flux-dependent definitions are recalled in Sections 2.3 and 2.4, respectively. In Section 2.5, growthdependent definitions are introduced. Section 2.6 is devoted to dead reactions, blocked reactions and their relationship. Section 2.7 shows how the integration of dynamic constraints affects the computation of vulnerabilities in Plasmodium falciparum model. Finally, Section 3 reviews the implemented tool CONTRABASS and Section 4 concludes the article.

Constraint-based models
This subsection recalls definitions (Oarga et al., 2020) and methods that will be used in the next sections.
Definition 2.1. A constraint-based model (Orth et al., 2011;Varma and Palsson, 1994) is a tuple fR; M; S; L; Ug where R is a set of reactions, M is a set of metabolites, S 2 R jMjÂjRj is the stoichiometric matrix, and L; U 2 R jRj are lower and upper flux bounds of the reactions.
Without loss of generality, it is assumed that L½r U½r 8r 2 R.
All reactions are associated with a set of reactant metabolites and a set of product metabolites (one of these two sets can be empty). For example, the reaction r : A ! 2B has a reactant metabolite A, and a product metabolite B with stoichiometric weight 2, i.e. two molecules of type B are produced per each molecule of type A that is consumed by r. The stoichiometric matrix S accounts for all the stoichiometric weights of the reactions, i.e. S½m; r is the stoichiometric weight of metabolite m 2 M for reaction r 2 R. Thus, • if S½m; r < 0 then m is a reactant and is consumed when r occurs. • if S½m; r > 0 then m is a product and is produced when r occurs. • if S½m; r ¼ 0 then m is neither consumed nor produced when r occurs. Constraint-based models can be expressed as bipartite graphs with two types of nodes representing reactions and metabolites. Hence, constraint-based models can be represented graphically as Petri nets (Heiner et al., 2008;Murata, 1989), where places, drawn as circles, model metabolites, and transitions, drawn as squares, model reactions. The presence of an arc from a place (transition) to a transition (place) means that the place is a reactant (product) of the reaction modelled by the transition. The weights of the arcs of the Petri net account for the stoichiometry of the constraint-based model. In other words, the stoichiometric matrix of a constraint-based model and the incidence matrix of its corresponding Petri net coincide.
Example 2.1.The Petri net in Figure 1 represents a simple constraintbased model that consists of 10 reactions and 7 metabolites. As an example, transition r 7 models the reaction r 7 : m e ! 2m f . Flux balance analysis (FBA) (Orth et al., 2010) is a mathematical procedure for the estimation of steady-state fluxes in constraintbased models. Among other possibilities, FBA can be used to predict the growth rate of an organism or the rate of production of a given metabolite.
Let v 2 R jRj be the vector of fluxes of reactions and v½r denote the flux of reaction r. The system is assumed to be mass balanced at steady state, this is, the rate of production and rate of consumption of metabolites is constant. This constraint is given by the expression S Á v ¼ 0, where S is the stoichiometric matrix. The steady-state mass balance fluxes of reactions, v, must also satisfy the lower and upper bounds L and U. Thus, the linear programming problem (LPP) for FBA is: where z 2 R jRj expresses the objective function. Let r g be the reaction that models growth (or biomass production) in a constraint-based model. Without loss of generality, it will be assumed that L½r g ! 0. An optimistic estimate for the growth rate of the modelled microorganism can be obtained by: The solution of the above LPP (Equation (2)) will be denoted by l max .

Essential reactions
A reaction is said to be essential if it is required by the organism to grow. In other words, the deletion of an essential reaction implies null growth. Consequently, these reactions have the potential to cause the death of the modelled organism.
Definition 2.2. (Oyelade et al., 2018) A reaction r 2 R is an essential reaction if the solution of the following LPP: max v½r g st: The set of essential reactions, which is denoted ER, can be computed straightforwardly by solving Equation (3) for each r 2 R.
Example 2.2. Let us say that in the Petri net of Figure 1 reaction r g models growth. In this model, reactions r 1 and r 2 are essential reactions. This is because, if the flux of any of these reactions is set to 0, then it is not possible to produce metabolite m b , which is necessary for the growth reaction r g .

Structural definitions
This section recalls how structural definitions of chokepoints and dead-end metabolites can be defined by making use of Petri net notation. Let X denotes a vertex of the net, i.e. a reaction or a metabolite. Then, • X(X • ) denotes the set of input (output) vertices of X. For instance, for a given reaction r 2 R; • r denotes its set of reactants and r • its set of products; for a given metabolite m 2 M; • m denotes its set of producing reactions and m • its set of consuming reactions.
A chokepoint is a reaction that is the only producer or the only consumer of a metabolite. More formally: Definition 2.3. (Chukualim et al., 2008;Yeh et al., 2004) The set of chokepoint reactions will be denoted as CP. Notice that the inhibition of a chokepoint reaction can lead either to the infinite accumulation of a substrate (which will not be used as expected or might be toxic) or to the depletion of a metabolite (which might be essential for the cell) (Chukualim et al., 2008;Yeh et al., 2004). The identification of chokepoints is, therefore, relevant to the identification of potential drug targets.
A dead-end metabolite is a metabolite without input or output reactions: Dead-end metabolites can indicate a potential shortcoming or incompleteness of the model since their concentration can only increase or decrease.

Flux-dependent definitions
The previous definitions only take into account the structure of the network and disregard the flux bounds of the reactions. In order to capture the fact that reactions can proceed forwards or backwards, e.g. reversible reactions, new sets of reactants, products, consumers and producers that take into account flux bounds were defined (Oarga et al., 2020): • Flux-dependent set of reactants of r: ? r ¼ fm 2 MjðSðm; rÞ < 0^U½r > 0Þ _ ðSðm; rÞ > 0^L½r < 0Þg • Flux-dependent set of products of r: r ? ¼ fm 2 MjðSðm; rÞ > 0^U½r > 0Þ _ ðSðm; rÞ < 0^L½r < 0Þg • Flux-dependent set of producers of m: ? m ¼ fr 2 Rjm 2 r ? g • Flux-dependent set of consumers of m: m ? ¼ fr 2 Rjm 2 ? rg Flux-dependent chokepoints and dead-end metabolites can be defined accordingly (Oarga et al., 2020): Example 2.3. In the Petri net in Figure 1, r 1 is a producer of m a , i.e. r 1 2 • m a ; r 2 is a flux-dependent chokepoint because it is the only consumer of m a , i.e. m a • ¼ fr 2 g and r 2 2 CP ? ; and m g is a flux-dependent dead-end metabolite, i.e. m g 2 DEM ? . In addition to chokepoints, flux bounds can be used to classify reactions as dead, reversible or non-reversible: Definition 2.9. A reaction r 2 R is non-reversible if r is not dead and r is not reversible.
From the above definitions, it can be deduced that r is nonreversible if ð0 L½r^0 < U½rÞ _ ðL½r < 0^U½r 0Þ. Notice that dead reactions will never have a non-zero flux. Since this might indicate a deficiency of the model or a blocked pathway of the organism, special attention will be paid to dead reactions.
The sets of dead, reversible and non-reversible reactions are denoted DR, RR and NR, respectively. By definition, DR, RR and NR represent a partition of the set of reactions R, i.e.
For the Petri net representation of a constraint-based model, dead reactions, reversible reactions and non-reversible reactions will be represented as rectangles with a cross inside, as double rectangles, and as rectangles, respectively.
Example 2.4. The Petri net in Figure 2a represents a constraint-based model with five metabolites and eight reactions. In this model, reactions r 3 , r 4 are reversible reactions (i.e. RR ¼ fr 3 ; r 4 g), reaction r 7 is a dead reaction (i.e. DR ¼ fr 7 g) and reactions r 1 ; r 2 ; r 5 ; r 6 ; r g are non-reversible reactions (i.e. NR ¼ fr 1 ; r 2 ; r 5 ; r 6 ; r g g).

Growth-dependent definitions
In contrast to the previous flux-dependent definitions of reactions, this section introduces growth-dependent definitions, i.e. the sets of reactions that depend on the growth constraints that are imposed on the model.

Essential reactions
Similar to essential reactions, which are those reactions that are necessary to produce non-null growth on the model, growthdependent essential reactions are those reactions that are necessary to produce a certain growth on the model. This certain growth will be expressed as c Á l max where c 2 ½0; 1 and l max is the solution of Equation (2).
A reaction is said to be a growth-dependent essential reaction for a given growth c Á l max if its deletion implies that the maximum possible growth is below c Á l max . More formally, Definition 2.10. Let l max be the solution of the LPP in Equation (2). Given c 2 ½0; 1, a reaction r 2 R is a growth-dependent essential reaction if the solution of LPP (Equation (3)) is lower than c Á l max or the LPP is infeasible.
The set of growth-dependent essential reactions for a given growth specified by c 2 ½0; 1 will be denoted ER c . This set can be computed straightforwardly by solving LPP (Equation (3)) for each reaction.
Special attention is given to the set of reactions ER 1 , as it will consist of those reactions that are necessary to produce the optimal growth of the model. This set will be named essential reactions for optimal growth (EROG).
Example 2.5. In the Petri net of Figure 1 where r g models growth, reactions r 1 ; r 2 ; r 3 ; r 4 are EROG (i.e. r 1 ; r 2 ; r 3 ; r 4 2 EROG). Reactions r 1 , r 2 are essential reactions, and thus, they are also EROG. If one, or both reactions, r 3 , r 4 are forced to have flux equal to 0, metabolite m, which is essential for growth, can only be produced through a non-optimal path. Then, the model will not be able to achieve optimal growth. Therefore, these reactions are EROG.

Dead, reversible and chokepoint reactions
Flux variability analysis (FVA) (Mahadevan and Schilling, 2003) is a mathematical procedure to compute the minimum and maximum fluxes of reactions that are compatible with some state. Let l max be the maximum growth calculated by FBA (Equation (2)). FVA is computed by solving two independent LPPs per reaction r 2 R. One programming problem maximizes the flux of r, v½r, and the other minimizes v½r. The constraints of both problems are the same: the steady state condition S Á v ¼ 0, the flux bounds L½r v½r U½r, and the maintenance of l max to a certain degree c 2 ½0; 1. The two programming problems for a given reaction r 2 R can be expressed as: where r g is the reaction that models growth. The computation of the flux bounds by means of FVA (Equation (4)), can be carried out in an optimal state, i.e. c ¼ 1, or in a suboptimal state, i.e. 0 c < 1. In the optimal state, all fluxes must be optimally directed towards growth, whereas in suboptimal states, fluxes are allowed to deviate towards other functionalities.
Let lb c ; ub c 2 R jRj be the result of computing FVA (Equation (4)) on a constraint-based model fR; M; S; L; Ug for a given c, i.e. lb c ½r and ub c ½r are the minimum and maximum fluxes given by FVA for reaction r. If the flux bounds L, U of the constrained-based model are replaced by lb c ; ub c , a new constraint-based model, fR; M; S; lb c ; ub c g, with more restrictive (and realistic) flux bounds is obtained.
Given c 2 ½0; 1, the sets of growth-dependent products, reactants, consumers, and producers of the model fR; M; S; L; Ug, which are denoted r c ; c r; m c ; c m, respectively, are defined as the fluxdependent products, reactants, consumers and products of fR; M; S; lb c ; ub c g as discussed in Section 2.4. Similarly, given fR; M; S; L; Ug and c 2 ½0; 1, the sets of growth-dependent dead, reversible and non-reversible reactions, which are denoted DR c ; RR c and NR c , the sets of growthdependent chokepoint reactions and dead-end metabolites, which are denoted as CP c and DEM c , are also defined as the corresponding flux-dependent elements of fR; M; S; lb c ; ub c g.
Example 2.6. Let us assume that r g in Figure 2a represents growth, i.e. the component of z in Equation (1) that corresponds to r g is equal to 1 and the rest of components of z are 0. In order to obtain growthdependent sets, the flux bounds computed by FVA with c ¼ 1 are assigned to the reactions and the net in Figure 2c is obtained (notice that other values of c can be considered). In this new net, r 2 , r 3 , r 5 , r 6 and r 7 are dead reactions, i.e. r 2 ; r 3 ; r 5 ; r 6 ; r 7 2 DR 1 , and r 4 , which was a reversible reaction in Figure 2a, becomes a non-reversible reaction, i.e. r 4 2 NR 1 .
In the model of Figure 2a, reaction r 3 is the only flux-dependent consumer of metabolite m b and, thus, is a flux-dependent chokepoint reaction, i.e. r 3 2 CP ? . However, in the growth-dependent model of Figure 2c, r 3 becomes a dead reaction, i.e. r 3 2 DR 1 , and hence, it is no longer considered a chokepoint reaction, i.e. r 3 6 2 CP 1 .
On the other hand, reaction r 4 is not the only producer or consumer of any metabolite in the model of Figure 2a, and thus, it is not a flux-dependent chokepoint reaction, i.e. r 4 6 2 CP ? . However, in the growth-dependent model of Figure 2c, reactions r 2 , r 3 and r 5 become dead reactions and reaction r 4 becomes the only consumer of m a and the only producer of m c , consequently, r 4 becomes a growthdependent chokepoint reaction r 4 2 CP 1 .
Given their importance in model tuning and drug discovery, the following sections focus on dead reactions, essential reactions and chokepoint reactions.

Dead reactions
This subsection explores the relationship between dead reactions and blocked reactions, and shows that the set of growth-dependent dead reactions is the same for any suboptimal state, i.e. for any growth strictly lower than the maximum growth l max . Such a set coincides with the set of blocked reactions. Recall that a reaction r 2 R is said to be blocked if its flux is 0 at any possible steady state. More formally: Definition 2.11. (Vlassis et al., 2014) A reaction r 2 R is a blocked reaction if for every v 2 R jRj such that S Á v ¼ 0, it holds v½r ¼ 0.
Example 2.7. In the Petri net in Figure 2a, reaction r 6 is a blocked reaction. This is because r 6 consumes m d , which is a DEM. If the flux of r 6 were positive (negative), the amount of m d would decrease(increase) indefinitely, which contradicts the steady-state constraint, therefore the only possible steady-state flux for r 6 is v½r 6 ¼ 0.
In Vlassis et al. (2014) blocked reactions are obtained by solving the following LPPs that compute the maximum and minimum feasible fluxes of the reactions subject to the steady state constraint A reaction with null maximum and minimum feasible flux is a blocked reaction. Note that this procedure is equivalent to computing FVA, see Equation (4), with c ¼ 0. Thus, the set of blocked reactions is equal to DR 0 .
Interestingly, the set of growth-dependent dead reactions, DR c is the same for any c such that 0 c < 1. In other words, the set of dead reactions in suboptimal states, regardless of the growth rate imposed on the model, is equivalent to the set of blocked reactions. This fact will be proved through several steps. Let us first prove that the range of feasible fluxes of a reaction r, i.e. the interval ½lb c ½r; ub c ½r, cannot increase as c increases.
Given c 2 ½0; 1, the range of feasible fluxes of r 2 R, ½lb c ½r; ub c ½r, is given by the solutions of FVA (Equation (4)). Notice that the constraints of such LPP define a convex set of possible solutions which can only shrink as c increases, i.e. as the constraint c Á l max z Á v becomes more restrictive. Thus, if c 1 < c 2 then lb c 1 ½r lb c 2 ½r and ub c 1 ½r ! ub c 2 ½r. h Then, the set of growth-dependent dead reactions cannot decrease with c: Proof. Let r 2 DR c 1 , i.e. lb c 1 ½r ¼ ub c 1 ½r ¼ 0, then by Lemma 2.1, it follows that lb c 2 ½r ¼ ub c 2 ½r ¼ 0 and hence r 2 DR c 2 . h Let us now show that the set of growth-dependent dead reactions cannot increase with c in suboptimal states, i.e. with c < 1: Proof. The convex set of possible solutions defined by the constraints in Equation (4) can only decrease as c increases, see Lemma 2.1. Assume there exist c a ; c b such that 0 c a < c b < 1; lb c a ½r < ub c a ½r and lb c b ½r ¼ ub c b ½r ¼ 0, i.e. r becomes dead when c a is increased to c b . Then, given that the constraints in Equation (4) are linear, if c b is further increased by 2 R such that > 0 and c b þ < 1, then lb c b þ ½r > ub c b þ ½r should hold which is not possible because a lower bound cannot exceed an upper bound. h From Lemmas 2.2 and 2.3, the following theorem can be derived straightforwardly: Theorem 2.4.DR c 1 ¼ DR c 2 8c 1 ; c 2 2 ½0; 1Þ.
Thus, in particular, the set of blocked reactions coincides with the set of dead reactions in suboptimal states.
Notice that DR 1 can be strictly greater than DR c with c 2 ½0; 1Þ. As an example, the next subsection analyses a constrained-based model in which DR 1 is strictly greater than DR c with c 2 ½0; 1Þ, see Figure 3d.

Case study: Plasmodium falciparum
This section presents the results obtained for the constraint-based model of P.falciparum (iAM-Pf480) (Abdel-Haleem et al., 2018). The model includes a total of 1083 reactions, 909 metabolites, and 480 genes. The sizes of the sets of flux-dependent reactions are jRRj ¼ 493; jNRj ¼ 590; jDRj ¼ 0. The number of initial fluxdependent chokepoint reactions is jCP ? j ¼ 453. The sets of essential reactions and EROG were also computed, and the following sizes were obtained: jERj ¼ 192 and jEROGj ¼ 317. Figure 3 shows the sizes of the growth-dependent sets NR c ; RR c ; CP c ; DR c ; ER c ; R À ER c À DR c in plots 3(a), 3(b), 3(c), 3(d), 3(e) and 3(f), respectively. To assess the impact of c on these sets, different values of c in the interval ½0; 1 have been used. In addition to the sizes of the sets obtained with c, the leftmost value (depicted in green) of plots 3(a), 3(b), 3(c) and 3(d) refers to the size of the flux-dependent set prior to FVA, i.e. NR, RR,CP ? and DR. Similarly, the leftmost value (in green) of plots 3(e) and 3(f) includes the sizes of the set of essential reactions ER and the set R À ER À DR, respectively.
Notice that if c ¼ 0, the constraint c Á l max z Á v in Equation (4) does not impose a minimum growth on the model, and only the steady state condition S Á v ¼ 0 must be satisfied. It can be seen that at c ¼ 0, the set of dead reactions exhibits an increase from the initial set of flux-dependent dead reactions: jDRj ¼ 0 and jDR 0 j ¼ 222. As noted in Section 2.6, this set DR 0 is equal to the set of blocked reactions, and keeps the same for every c 2 ½0; 1Þ, see Figure 3d. The second increase in the set of dead reactions takes place at c ¼ 1 and it is due to the fact that in the optimal growth state, the flux must be necessarily distributed through optimal paths for biomass production and no flux can be diverted through other paths. Thus, nonoptimal reactions for biomass production become dead reactions.
With respect to reversible reactions, the steady-state constraint S Á v ¼ 0 reduces the size of this set from jRRj ¼ 493 to jRR 0 j ¼ 210, see Figure 3b. Such a reduction is caused by the blocked reactions that belonged to RR and become dead reactions, and by the reversible reactions that become non-reversible with the steady-state constraint.
Similarly, blocked reactions that belonged to the set of fluxdependent non-reversible reactions NR become dead reactions with the steady-state constraint. However, due to the significant amount of reversible reactions that become non-reversible reactions, the size of the set increases from jNRj ¼ 590 to jNR 0 j ¼ 651, see Figure 3a.
With regard to chokepoint reactions, see Figure 3c, the number of flux-dependent chokepoints is jCP ? j ¼ 453. At c ¼ 0, there is a decrease to jCP 0 j ¼ 450, and then the set increases slowly until jCP 0:99 j ¼ 507. As in the sets of non-reversible reactions and reversible reactions, the set of chokepoints decreases at c ¼ 1 as many reactions become dead reactions.
Notice that the set of chokepoints at c ¼ 1, CP 1 , is smaller than the set of flux-dependent chokepoints, CP ? . Moreover, CP 1 is not contained in CP ? . This is due to the changes produced in the sets of non-reversible reactions and reversible reactions as c increases. Recall that a reversible reaction is considered simultaneously both as a consumer and as a producer of a given metabolite. If this reaction becomes non-reversible, then it might not be a producer and consumer of a given metabolite. This fact can lead to the appearance of new growth-dependent chokepoints in CP 1 that were not chokepoints in CP ? . On the other hand, dead reactions are neither consumers nor producers and, hence, they are not chokepoints. Thus, the increase of dead reactions in DR 1 can imply changes in CP 1 with respect to CP ? . Figure 3e reports the size of the set of essential reactions ER c with respect to c. Recall that the leftmost value of this graph refers to the ER set, i.e. the set of essential reactions with no growth constraints.
Notice that no reaction is mandatory to produce a null growth, hence at c ¼ 0, jER 0 j is always 0. Furthermore, for positive values of c the size of ER c increases as the size of RR c decreases. This is due to the fact that, as c increases, some reversible reactions become nonreversible, and the flux of these reactions is necessary, i.e. essential for growth. On the other hand, as expected, the amount of growthdependent essential reactions increases with c.
Finally, Figure 3f shows the size of the set R À ER c À DR c with respect to c, with the leftmost point corresponding to R À ER À DR. This set is of interest as it is composed of those reactions that can contribute to growth, that is, they are not dead reactions, and at the same time are not growth-dependent essential reactions, that is, if one of these reactions is knocked out the growth can still be kept the same (notice that this does not mean that all the reactions in this set can be knocked simultaneously without affecting growth). The fact that knocking out one of these reactions does not reduce the growth might then imply the presence of redundancy in the metabolism. Therefore, this set of reactions provides the metabolism with resilience and flexibility. As it can be seen, the size of this set decreases with c. This fact means that producing higher growths is more demanding as it requires a higher amount of reactions to sustain it. In other words, the higher the growth the more reactions are necessary to produce such growth, and therefore, less flexibility is given to the metabolism.

Implementation
CONTRABASS is a software tool for the computation of vulnerabilities in GEMs. CONTRABASS is distributed as a Python command line tool but can also be executed through an online web server at http://contrabass.unizar.es. The CONTRABASS web server is designed to offer an intuitive interface to access the operations of the tool.
The tool takes as an input a model in systems biology markup language (SBML) (Hucka et al., 2019) of Levels 1, 2 or 3, and computes the set of chokepoints reactions, essential reactions, dead reactions and dead-end metabolites on the model by taking into account the dynamic constraints for the model as explained in the article. The results are then exported as a spreadsheet file and as an interactive HTML report. The operations that CONTRABASS allows include the computation of sets of chokepoints, dead, reversible, nonreversible and essential reactions with different values of c; computation and removal of dead-end metabolites from a model; and the update of the flux bounds of the reactions according to FVA.
In addition to the above, through the interactive HTML report users can also access the data available in the model, this is, reactions, genes and metabolites along with their databases identifiers if (e) ( f) Fig. 3. Sizes of the sets of reactions NRc; RRc; CPc; DRc; ERc and R À ERc À DRc of P.falciparum for c 2 ½0; 1. The leftmost value of each plot corresponds to NR, RR, CP ? , DR, ER and R À ER À DR, respectively available; explore the reaction sets defined in this article; and also explore the intersection of different sets of critical reactions. The documentation of the tool is available at https://contrabass. readthedocs.io. CONTRABASS makes use of the Python toolbox COBRApy (Ebrahim et al., 2013). The source code of the command line tool and the web server is available at https://github.com/ openCONTRABASS/CONTRABASS and https://github.com/ openCONTRABASS, respectively. All the code is released under GPL-3.0 License.

Discussion
The development of new drugs is an expensive and costly process. This process needs to be streamlined in order to face current health threats as those posed by MDR. In this work, computational methods to identify critical reactions, i.e. vulnerabilities, in the metabolism of microorganisms are proposed. When applied on models of pathogenic bacteria, such critical reactions are potential drug targets.
The developed methods classify the reactions of the models in different sets, e.g. essential reactions, chokepoint reactions, blocked reactions, etc., according to their role in the metabolism. In order to exploit the information available in the model as much as possible, flux-and growth-dependent sets of reactions have been defined. Such sets capture the role of the reactions as a function of the growth rate of the organism. Thus, it is possible to unveil how the cell uses its reactions to achieve its growth. The relationship among these sets has been explored and it has been shown that the set of dead reactions keeps constant for any non-optimal growth rate. Moreover, the computation of other flux-and growth-dependent sets, such as dead-end metabolites and dead reactions, can help in the identification of potential shortcomings of the model. In order to facilitate their use, an online web server and a command line tool implementing the discussed methods have been made publicly available.